CompDay 11 - Colin Fisher

#Importing packages

#Numpy for vector calculations and plotting
import numpy as np

#Matplotlib for plotting (and animating)
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

#Sympy for vector calculations
import sympy as sym
sym.init_printing()

#Lambdify for converting sympy functions/values to numpy ones
from sympy.utilities.lambdify import lambdify

Part A

Written Work:

After determining the sum of forces acting on each mass (as seen above), we get the following matrices for $\textbf{K}$ and $\textbf{M}$:

#Defining matrices K and M

#Defining symbols
k, m = sym.symbols('k, m')

#Defining matrix M
M = sym.Matrix([[m, 0, 0, 0],[0, m, 0, 0],[0, 0, m, 0],[0, 0, 0, m]])

M #displaying result

$\displaystyle \left[\begin{matrix}m & 0 & 0 & 0\\0 & m & 0 & 0\\0 & 0 & m & 0\\0 & 0 & 0 & m\end{matrix}\right]$

The above result is our $\textbf{M}$ matrix for our system

#Defining matrix K
K = sym.Matrix([[2*k, -k, 0, 0],[-k, 2*k, -k, 0],[0, -k, 2*k, -k],[0, 0, -k, 2*k]])

K #displaying result

$\displaystyle \left[\begin{matrix}2 k & - k & 0 & 0\\- k & 2 k & - k & 0\\0 & - k & 2 k & - k\\0 & 0 & - k & 2 k\end{matrix}\right]$

The above result is our $\textbf{K}$ matrix for our system

Part B

We can then calculate the normal frequencies through taking the determinant of $\textbf{K}-\omega^2\textbf{M}$, setting the result equal to zero, and solving for $\omega^2$:

#Calculating normal frequencies

#Defining new symbol
ω = sym.symbols('ω')

#Defining our relationship
KwM = K - ω**2*M
KwM #Checking result

$\displaystyle \left[\begin{matrix}2 k - m ω^{2} & - k & 0 & 0\\- k & 2 k - m ω^{2} & - k & 0\\0 & - k & 2 k - m ω^{2} & - k\\0 & 0 & - k & 2 k - m ω^{2}\end{matrix}\right]$

#Calculating determinant of above matrix
det = sym.det(KwM)
total_norm_freq = sym.solve(det, ω**2) #Solving for ω^2

#Saving these determinants as variables
ω_norm1 = total_norm_freq[0]
ω_norm2 = total_norm_freq[1]
ω_norm3 = total_norm_freq[2]
ω_norm4 = total_norm_freq[3]

total_norm_freq #Displaying result

$\displaystyle \left[ \frac{k \left(\frac{3}{2} - \frac{\sqrt{5}}{2}\right)}{m}, \ \frac{k \left(\frac{5}{2} - \frac{\sqrt{5}}{2}\right)}{m}, \ \frac{k \left(\frac{\sqrt{5}}{2} + \frac{3}{2}\right)}{m}, \ \frac{k \left(\frac{\sqrt{5}}{2} + \frac{5}{2}\right)}{m}\right]$

Thus our normal frequencies are: $$\omega^2_1 = \frac{k(3/2-\sqrt{5}/2)}{m}$$ $$\omega^2_2 = \frac{k(5/2-\sqrt{5}/2)}{m}$$ $$\omega^2_3 = \frac{k(5/2+\sqrt{5}/2)}{m}$$ $$\omega^2_4 = \frac{k(3/2+\sqrt{5}/2)}{m}$$

Part C

Plugging in each of these normal frequencies to the original equation $\textbf{K}-\omega^2\textbf{M}$, as outlined in Part b, we can calculate their corresponding eigenvectors to gather their corresponding normal mode:

#Defining new symbols
a1, a2, a3, a4 = sym.symbols('a_1, a_2, a_3, a_4')

#Defining general matrix - used in determining
A_matrix = sym.Matrix([a1,a2,a3,a4])

#First Normal Mode - substituting first solution from Part B into matrix, multiplying by general matrix, 
#and getting relationships between variables a1, a2, a3, and a4 based on resulting system of equations

arg1 = sym.solve(KwM.subs(ω**2,ω_norm1)*A_matrix,(a1,a2,a3,a4))
mode_norm1 = A_matrix.subs(arg1).subs(a4,1) #substituting above arguments into general matrix and evaluating when a4 = 1

mode_norm1 #Displaying result of first normal mode

$\displaystyle \left[\begin{matrix}1\\\frac{1}{2} + \frac{\sqrt{5}}{2}\\\frac{1}{2} + \frac{\sqrt{5}}{2}\\1\end{matrix}\right]$

#Second Normal Mode - same process as above
arg2 = sym.solve(KwM.subs(ω**2,ω_norm2)*A_matrix,(a1,a2,a3,a4)) #same method of generating arguments
mode_norm2 = A_matrix.subs(arg2).subs(a4,1) #same substitution

#Third Normal Mode - same process as above
arg3 = sym.solve(KwM.subs(ω**2,ω_norm3)*A_matrix,(a1,a2,a3,a4)) #same method of generating arguments
mode_norm3 = A_matrix.subs(arg3).subs(a4,1) #same substitution

#Fourth Normal Mode - same process as above
arg4 = sym.solve(KwM.subs(ω**2,ω_norm4)*A_matrix,(a1,a2,a3,a4)) #same method of generating arguments
mode_norm4 = A_matrix.subs(arg4).subs(a4,1) #same substitution

mode_norm1, mode_norm2, mode_norm3, mode_norm4 #Displaying ALL results of all normal modes

$\displaystyle \left( \left[\begin{matrix}1\\\frac{1}{2} + \frac{\sqrt{5}}{2}\\\frac{1}{2} + \frac{\sqrt{5}}{2}\\1\end{matrix}\right], \ \left[\begin{matrix}-1\\\frac{1}{2} - \frac{\sqrt{5}}{2}\\- \frac{1}{2} + \frac{\sqrt{5}}{2}\\1\end{matrix}\right], \ \left[\begin{matrix}1\\\frac{1}{2} - \frac{\sqrt{5}}{2}\\\frac{1}{2} - \frac{\sqrt{5}}{2}\\1\end{matrix}\right], \ \left[\begin{matrix}-1\\\frac{1}{2} + \frac{\sqrt{5}}{2}\\- \frac{\sqrt{5}}{2} - \frac{1}{2}\\1\end{matrix}\right]\right)$

The above result contains the normal modes for our system

Part D

#Animating 4 different normal modes

#Defining plot
fig, ax = plt.subplots(figsize=(10,10))

#Defining arrays for initial x positions AND normal modes from part c
x_init = np.array([-9, -3, 3, 9])
norm1 = np.array([mode_norm1[0], mode_norm1[1], mode_norm1[2], mode_norm1[3]])
norm2 = np.array([mode_norm2[0], mode_norm2[1], mode_norm2[2], mode_norm2[3]])
norm3 = np.array([mode_norm3[0], mode_norm3[1], mode_norm3[2], mode_norm3[3]])
norm4 = np.array([mode_norm4[0], mode_norm4[1], mode_norm4[2], mode_norm4[3]])

#Converting sympy functions to numpy functions of normal frequencies from part b
ω1_func = lambdify((k, m), ω_norm1, 'math')
ω1 = ω1_func(1,1) #plugging in k = 1 and m = 1 to solve for ω1
ω2_func = lambdify((k, m), ω_norm2, 'math')
ω2 = ω2_func(1,1) #plugging in k = 1 and m = 1 to solve for ω2
ω3_func = lambdify((k, m), ω_norm3, 'math')
ω3 = ω3_func(1,1) #plugging in k = 1 and m = 1 to solve for ω3
ω4_func = lambdify((k, m), ω_norm4, 'math')
ω4 = ω4_func(1,1) #plugging in k = 1 and m = 1 to solve for ω4

#Plotting each mass
m1_norm1, = ax.plot(x_init[0],1,'s',color='red', markersize = '25') #Plot a single point representing m1 of normal mode 1
m2_norm1, = ax.plot(x_init[1],1,'s',color='red', markersize = '25') #Plot a single point representing m2 of normal mode 1
m3_norm1, = ax.plot(x_init[2],1,'s',color='red', markersize = '25') #Plot a single point representing m2 of normal mode 1
m4_norm1, = ax.plot(x_init[3],1,'s',color='red', markersize = '25') #Plot a single point representing m3 of normal mode 1

m1_norm2, = ax.plot(x_init[0],2,'s',color='blue', markersize = '25') #Plot a single point representing m1 of normal mode 2
m2_norm2, = ax.plot(x_init[1],2,'s',color='blue', markersize = '25') #Plot a single point representing m2 of normal mode 2
m3_norm2, = ax.plot(x_init[2],2,'s',color='blue', markersize = '25') #Plot a single point representing m2 of normal mode 2
m4_norm2, = ax.plot(x_init[3],2,'s',color='blue', markersize = '25') #Plot a single point representing m3 of normal mode 2

m1_norm3, = ax.plot(x_init[0],3,'s',color='green', markersize = '25') #Plot a single point representing m1 of normal mode 3
m2_norm3, = ax.plot(x_init[1],3,'s',color='green', markersize = '25') #Plot a single point representing m2 of normal mode 3
m3_norm3, = ax.plot(x_init[2],3,'s',color='green', markersize = '25') #Plot a single point representing m2 of normal mode 3
m4_norm3, = ax.plot(x_init[3],3,'s',color='green', markersize = '25') #Plot a single point representing m3 of normal mode 3

m1_norm4, = ax.plot(x_init[0],4,'s',color='purple', markersize = '25') #Plot a single point representing m1 of normal mode 4
m2_norm4, = ax.plot(x_init[1],4,'s',color='purple', markersize = '25') #Plot a single point representing m2 of normal mode 4
m3_norm4, = ax.plot(x_init[2],4,'s',color='purple', markersize = '25') #Plot a single point representing m2 of normal mode 4
m4_norm4, = ax.plot(x_init[3],4,'s',color='purple', markersize = '25') #Plot a single point representing m3 of normal mode 4

#Defining animation function for time t
def animate(t):
    #Calculating positions
    x_norm1 = x_init + norm1*(np.cos(ω1*t)+np.sin(ω1*t)) #positions for all masses of normal mode 1
    x_norm2 = x_init + norm2*(np.cos(ω2*t)+np.sin(ω2*t)) #positions for all masses of normal mode 2
    x_norm3 = x_init + norm3*(np.cos(ω3*t)+np.sin(ω3*t)) #positions for all masses of normal mode 3
    x_norm4 = x_init + norm4*(np.cos(ω4*t)+np.sin(ω4*t)) #positions for all masses of normal mode 4
    
    #Updating positions
    #Normal Mode 1
    m1_norm1.set_xdata(x_norm1[0])
    m2_norm1.set_xdata(x_norm1[1])
    m3_norm1.set_xdata(x_norm1[2])
    m4_norm1.set_xdata(x_norm1[3])
    #Normal Mode 2
    m1_norm2.set_xdata(x_norm2[0])
    m2_norm2.set_xdata(x_norm2[1])
    m3_norm2.set_xdata(x_norm2[2])
    m4_norm2.set_xdata(x_norm2[3])
    #Normal Mode 3
    m1_norm3.set_xdata(x_norm3[0])
    m2_norm3.set_xdata(x_norm3[1])
    m3_norm3.set_xdata(x_norm3[2])
    m4_norm3.set_xdata(x_norm3[3])
    #Normal Mode 4
    m1_norm4.set_xdata(x_norm4[0])
    m2_norm4.set_xdata(x_norm4[1])
    m3_norm4.set_xdata(x_norm4[2])
    m4_norm4.set_xdata(x_norm4[3])
    
    #Returning the motion of all 16 masses
    return m1_norm1, m2_norm1, m3_norm1, m4_norm1, m1_norm2, m2_norm2, m3_norm2, m4_norm2, m1_norm3, m2_norm3, m3_norm3, m4_norm3, m1_norm4, m2_norm4, m3_norm4, m4_norm4

#Creating list of time values in seconds
tlist = np.linspace(0,20,200)

#Running animation function
ani = FuncAnimation(fig, animate, tlist, interval=50, blit=True) 

#Plotting and animating!
plt.title('Four Normal Oscillation Modes of Four Equal Masses Connected by Equal Springs')
plt.xlabel('X Position (m)')
plt.xlim(-12,12)
plt.ylim(0,5)

#Giving new tick marker values for y-axis - code modified from https://stackoverflow.com/questions/11244514/modify-tick-label-text
fig.canvas.draw()
labels = [item.get_text() for item in ax.get_yticklabels()]
labels[0] = ''
labels[1] = 'Normal Mode 1'
labels[2] = 'Normal Mode 2'
labels[3] = 'Normal Mode 3'
labels[4] = 'Normal Mode 4'
labels[5] = ''
ax.set_yticklabels(labels)

#Closing plot and converting HTML
plt.close()
HTML(ani.to_jshtml())

Once Loop Reflect